function estimate_pf_CD(treatment,k)  

%%% Settings
nMinObs = 25; % minimum number of obs required for GMM


% backward-looking rolling window for data (no. of quarters incl. current)
movingWindowSize = 20; 

% optimizer options
optionsUnc = optimoptions(@fminunc,'Display','off','StepTolerance',1e-14,'FunctionTolerance',1e-14); 
optionsNM = optimset('Display','off'); 

%%% load data
dataIn = readtable(['../../B_temp/pf_input_naicsq' num2str(k) '_' treatment '.csv']);

% create lists of industries and dates  
allNaicsCodes   = unique(dataIn{:,'naicsq'});   
allDates        = unique(dataIn{:,'date'});   
 
%%% estimation
 
% loop over industry-quarters ... 
T = length(allDates);
N = length(allNaicsCodes); 

params = [.85;.15;0];

for i = 1:N
    
    disp(['industry ' num2str(allNaicsCodes(i)) ' ... '])

    for t = 1:T  
        
        idx = (i-1)*T + t;
        paramsPre = params;
        
        % industry-quarter considered
        thisNaicsCode = allNaicsCodes(i);
        thisDate = allDates(t);
        
        % save estimates: indices
        naicsList(idx,1) = thisNaicsCode;
        datesList{idx,1} = thisDate;

        % select data 
        relevantRows = (dataIn.naicsq == thisNaicsCode & thisDate - dataIn.date >= 0 & thisDate - dataIn.date <= movingWindowSize);  
        
        % start estimation if observations available   
        if sum(relevantRows) < nMinObs 
            
            thetaCogsList(idx,1) = NaN;
            thetaCapitalList(idx,1) = NaN; 
            nObsList(idx,1) = sum(relevantRows);
            
        else
            
            % data  
            data = dataIn{relevantRows,{'y','v','k'}};  
                
            % first-stage measurement error correction
            dataPoly = ones(size(data(:,2))); 
            for ii=1:3
                dataPoly = [dataPoly data(:,2).^ii];
                dataPoly = [dataPoly data(:,3).^ii];
                for jj=1:3
                    dataPoly = [dataPoly data(:,2).^ii .* data(:,3).^jj];
                end
            end     
                
            coeff = dataPoly\data(:,1); 
            y_corr = dataPoly*coeff;
            Ly_corr = [NaN; y_corr(1:end-1)];
            
            % create data with lags in panel structure
            id = dataIn{relevantRows,{'id'}}; 
            tt = dataIn{relevantRows,{'date'}}; 
 
            data = [y_corr, ...
                    dataIn{relevantRows,{'v','k'}}, ...
                    Ly_corr,...
                    dataIn{relevantRows,{'Lv','Lk'}},...  
                    dataIn{relevantRows,{'L4v'}},... 
                    id,...
                    [NaN; id(1:end-1)],...
                    tt,...
                    [NaN; tt(1:end-1)],...
                    ]; 
                
            data = data(data(:,end-1)==data(:,end)+1 & data(:,end-3) == data(:,end-2),1:end-4);

            % again require minimal number of observations
            if size(data,1) < nMinObs 
                naicsList(idx,1) = thisNaicsCode;
                datesList{idx,1} = thisDate;
                thetaCogsList(idx,1) = NaN;
                thetaCapitalList(idx,1) = NaN; 
            else
            
            % OLS 
            paramsStartOLS = ols(data(:,1),[data(:,2), data(:,3)]);
            paramsStart = paramsStartOLS(1:2); 
            
            % GMM
                
            allFvals = []; 

            % unconstrained
            [paramsUnc1, fvalUnc1, flagUnc1] = fminunc(@(params) GMM_objective_CD(params, data), paramsStartOLS, optionsUnc);
            allFvals = [allFvals  fvalUnc1]; 

            [paramsUnc2, fvalUnc2, flagUnc2] = fminunc(@(params) GMM_objective_CD(params, data), paramsPre, optionsUnc);
            allFvals = [allFvals  fvalUnc2];  

            % nelder-mead    
            [paramsNM, fvalNM, flagNM] = fminsearch(@(params) GMM_objective_CD(params, data), paramsStartOLS, optionsNM);
            allFvals = [allFvals fvalNM];
            
            % take the best estimate over all minimizers
            kk=0;
            for minimizer = {'Unc1','Unc2','NM'}
                kk=kk+1;
                if eval(['fval' minimizer{1}]) <= min(allFvals)
                    params = eval(['params' minimizer{1}]);
                    fval   = eval(['fval' minimizer{1}]);
                    flag   = eval(['flag' minimizer{1}]);
                    method = minimizer{1};
                end
            end
            
            thetaCogsList(idx,1) = params(1); 
            thetaCapitalList(idx,1) = params(2);   
          
    end
    end
    end 
end

%% results
results = table(naicsList,datesList,thetaCogsList);
results.Properties.VariableNames = {'naicsq', 'date_quarterly', 'thetaCogsCD'}; 
writetable(results,['../../B_temp/pf_CD_naicsq' num2str(k) '_'  treatment '.csv'],'Delimiter',',') 

end